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1.  General 


'y  The  long-term  goal  of  our  research  at  the  Weizmann  Institute  is  the  devel¬ 
opment  of  multi-level  methods  for  solving  all  types  of  large-scale  problems  in  sci¬ 
ence  and  engineering.  In  addition  to  an  extensive  and  diverse  development  of 
multigrid  solvers  for  differential  and  integral  equations,  completely  new  multilevel 
approaches  have  recently  been  introduced  to  the  areas  of  large-scale  global  opti¬ 
mization,  statistical  physics  and  calculation  of  many-body  interactions  (see  our 
review,[l]). 

The  three-year  contract  summarized  below  has  undertaken  two  multigrid 
themes:  grid  adaptation  and  time-dependent  problems.  Due  to  some  personnel 
reasons,  it  turned  out  that  during  the  first  year  and  part  of  the  second  the  main 
advance  had  been  in  the  first  theme,  while  in  the  second  phase  effort  has  been 
concentrated  in  the  latter.  The  paper  summarizing  our  work  on  grid  adaptation 
has  already  appeared  [2],  while  the  final  paper  on  time-dependent  problems  is  still 
in  preparation,  and  will  be  sent  later.  For  scientific  reports  attached  herein  — 
see  the  list  of  references. 


2.  Grid  adaptation 


r. 


r\<» 
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Since  our  research  on  multigrid  mesh  refinement  techniques  is  fully  reported 
in  [2],  we  confine  the  present  description  to  a  brief  summary. 

The  purpose  of  our  study  was  to  provide  criteria  for  optimizing  meshsizes  near 
singularities  and  to  develop  fast  and  flexible  multigrid  methods  for  creating  the 
nonuniform  grids,  their  difference  equations  and  their  solutions.  For  simplicity,  the 
Poisson  problem  was  studied,  with  singularities  introduced  either  in  the  forcing 
terms  (algebraic  singularities  or  sources)  or  in  the  shape  of  the  boundaries  (re¬ 
entrant  comers).  Local  refinements  were  created  by  multigrid  structures  in  which 
some  extra  finer  leyfels  cover  increasingly  narrower  neighborhoods  of  the  singularity, 
as  proposed  in-{3f.  The  main  innovations  developed  under  the  present  contract 
are  described  below. 


2.1  Local  relaxation  sweeps 

Structural  singularities,  such  as  reentrant  comers  or  discontinuities  in  co¬ 
efficients  (of  the  interior  or  boundary  differential  operators)  result  in  degraded 
multigrid  performance  for  the  following  two  reasons:  1)  Error  components  slow 
to  converge  in  relaxation  are  approximate  solutions  of  the  homogeneous  differ¬ 
ential  equation.  Near  a  structural  singularity  components  of  this  type  can  have 
large  derivatives  and  hence  large  errors  in  their  interpolation  from  any  coarse- grid 


approximation.  2)  When  transferring  residuals  from  a  fine  grid  to  a  coarse  one, 
unusual  weighting  should  be  given  to  residuals  near  the  singularity. 

Experiments  show  that  the  asymptotic  convergence  degrades  only  when  V 
cycles  (not  when  W  cycles)  axe  used,  and  especially  when  the  number  of  levels  is 
large.  Which  is  why  this  trouble  particularly  bothers  us  when  local  refinements 
(needed  near  structural  singularities  for  accuracy  purposes)  are  used:  With  local 
refinements  many  levels  are  usually  introduced,  and  W  cycles  are  prohibitively 
expensive  (since  the  number  of  points  on  coarse  levels  is  not  necessarily  small 
compared  to  that  on  fine  levels). 

We  have  therefore  developed  a  general  approach  to  obtain  the  full  efficiency 
of  multigrid  V  cycles  despite  structural  singularities.  The  technique  is  to  make 
special  local  relaxation  passes  over  a  small  number  of  points  near  the  singularity. 
This  eliminates  the  large  interpolation  errors  (which  are  large  only  locally  there) 
and  reduces  large  local  residuals,  hence  making  their  correct  weighting  immaterial. 

Experiments  with  this  approach  were  conducted  on  the  Poisson  problem  with 
Dirichlet  boundary  conditions  in  a  domain  with  a  2rr  re-entrant  boundary.  They 
showed  that  the  same  asymptotic  convergence  rates  as  in  regular  problems  (i.e.,  as 
in  the  Poisson  problem  on  regular  domains)  can  be  obtained  by  applying,  on  each 
grid  at  its  turn,  one  local  relaxation  pass  over  a  few  points  near  the  singularity 
before  each  usual  relaxation  sweep  over  the  entire  grid.  The  number  of  points 
where  the  local  relaxation  should  be  performed  is  a  small  fraction  of  the  entire 
grid,  hence  the  amount  of  extra  computations  involved  in  the  special  passes  is 
negligible  compared  to  the  usual  multigrid  work. 

Incidentally,  the  technique  of  extra  relaxation  sweeps  over  a  small  region  near 
the  boundary  was  later  incorporated  into  a  general  rigorous  theory  of  multigrid 
solvers.  It  is  used  even  in  some  cases  without  structural  singularities,  to  ensure 
that  the  “interior  convergence  factor”  is  always  attained  (see  [4]). 


2.2  The  A-FMG  algorithm 

A  general  technique  for  treating  singularities  is  by  local  mesh  refinement.  As 
mentioned  above,  a  method  for  creating  local  refinements  as  a  natural  part  of 
multigrid  fast  solvers  had  been  suggested  in  [3,  §8]  and  [5,  §9],  and  was  further 
studied  by  us  as  part  of  this  contract. 

The  main  question  we  first  treated  is  how  to  modify  the  FMG  algorithm 
[5,  §7]  for  the  situation  where  many  levels  of  local  refinements  are  created.  In 
that  situation,  when  finer  levels  cover  much  smaller  subdomains,  the  number  of 
gridpoints  on  coarser  levels  is  not  necessarily  small  in  comparison  to  their  number 
on  finer  levels,  hence  the  work  expended  on  coarser  levels  by  the  usual  FMG 
algorithm  is  much  too  large. 


A  modified  algorithm,  called  A-FMG,  has  been  devised.  It  is  based  on  the 
same  exchange  rate  (the  Lagrange  multiplier  A)  used  in  the  grid-adaptation  equa¬ 
tions  (see  [3,  §8]).  For  each  value  of  A  these  equations  determine  a  set  of  local- 
refinement  grids  on  which  the  problem  should  be  discretized.  The  A-FMG  algo¬ 
rithm  uses  a  sequence  Ao  >  Aj  >  A2  >  . . .  of  exchange-rate  values,  where  Ao  is 
large  so  that  its  set  of  local-refinement  grids  is  relatively  coarse,  giving  a  crude 
solution  for  a  small  amount  of  work,  and  A,-  =  Q(A,-_i,  with  constant  factor  a  <  1. 
For  each  A,-  (*  >  1),  the  first  approximation  is  obtained  by  interpolating  the  A,_i 
solution  from  the  A,-_i  set  of  grids  into  the  At-  set.  Assuming  At_j  to  be  close 
enough  to  A,-,  this  first  approximation  is  such  that  a  multigrid  cycle  with  only  a 
small  number  of  relaxation  sweeps  on  each  grid  of  Aj  is  required  to  obtain  the  Aj 
solution  to  within  truncation  errors. 

A  detailed  analysis  shows  that  the  efficiency  of  the  algorithms  is  not  sensitive 
to  the  exact  value  of  a  =  Aj/A ,-_i.  For  example,  for  second-order  problems  in  two 
dimensions,  any  value  of  a  roughly  in  the  range  .02  <  a  <  .25  will  give  results 
close  to  optimal.  In  our  numerical  experiments  a  =  1/16  =  .0625  was  taken. 

The  numerical  experiments  were  first  carried  out  for  a  problem  with  right- 
hand  side  singularity,  in  order  to  separate  the  above  question  (the  efficiency  of 
FMG)  from  other  questions  (specific  to  structural  singularities;  see  §2.1).  Namely, 
we  have  solved  the  equation  Au  =  72r7~2  on  the  square  (0  <  x  <  Rq,  0  < 
V  <  -Ro}  where  r  is  the  distance  from  the  comer  (0,0).  The  Dirichlet  boundary 
conditions  were  such  that  the  exact  solution  is  tt  =  r7.  The  known  solution  allows 
us  to  make  precise  comparisons  of  algebraic  errors  with  truncation  errors. 

We  then  experimented  also  with  the  problem  with  2 ir  reentrant  comer,  aided 
by  local  relaxation  sweeps  as  described  above  (§2.1).  The  boundary  conditions 
were  chosen  so  that  the  exact  solution  of  the  differential  equations  was  u  = 
r1/2  sin</5/2,  where  (r,  <ys)  are  radial  coordinates  such  that  the  re-entrant  comer 
is  at  r  =  0  and  the  re-entrant  boundary  is  at  ip  =  0  and  <p  =  2x.  This  solution 
exhibits  the  typical  singularity  at  such  a  comer. 

The  grid-adaptation  equations  [3,  §8]  yield  for  each  of  these  problems  a  set 
of  grids  where  the  finest  meshsize  h  at  distance  r  from  the  singularity  is  given  by 
a  known  rule  h  =  h(r,  A).  For  example,  for  the  right-hand  singularity  problem 
h  =  0(A1/4r^2-7)/4).  We  have  used  for  each  A  a  set  of  grids  approximately 
satisfying  this  rule,  and  applied  the  A-FMG  algorithm  described  above.  Two 
important  results  were  established: 

1.  The  A-FMG  algorithm  solves  the  equations  in  linear  times;  i.e.,  for  each  prob¬ 
lem,  the  solution  time  is  proportional  to  the  total  number  of  gridpoints  on  all 
levels.  (The  usual  FMG  algorithm  falls  far  behind  such  linear  dependence  in 
the  case  under  study,  because  of  the  many  levels  of  local  refinements). 

2.  Due  to  the  optimal  choice  of  h  as  a  function  of  r  and  the  above  efficiency 
of  the  A-FMG,  each  problem  was  solved  in  the  same  over- all  accuracy-to- 
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work  convergence  rate  as  a  regular  problem.  Namely,  despite  the  singularity, 
the  error  E  in  approximating  the  true  differential  solution  decreases  as  IV-1, 
where  W  is  the  over-all  computational  work;  this  rate  E  =  0(W-1)  is  the  best 
possible  for  any  second-order  approximation  to  any  two-dimensional  problem. 


2.3  Interior  singularity.  Conservative  local  refinement 

Source  singularities  and  other  interior  singularities  appear  in  many  applica¬ 
tions.  As  a  model  we  took  the  equation  Au(i,y)  =  4nS,  where  6  is  the  Dirac’s 
distribution,  with  boundary  conditions  that  yield  the  solution  u  =  log(x2  +  y2). 
The  singularity  is  at  (0,0),  which  was  placed  at  the  center  of  the  square  domain. 
We  tested  local  refinement  techniques  on  this  model.  This  led  to  some  modifica¬ 
tions  of  those  techniques. 

First  we  found  that  no  local  refinement  is  really  needed  to  obtain  regular  ac¬ 
curacy  far  from  the  source:  if  the  second-order  differencing  is  in  conservation  form, 
and  if  the  correct  source  strength  is  represented,  then  solution  errors  are  always 
0(/i2/r2)  at  distance  r  from  the  source.  Furthermore,  adding  non-conservative 
local  refinements  very  substantially  increases  these  errors  —  proportionately  in 
fact  to  the  spurious  (numerically  created)  fluxes.  It  is  therefore  important  to  use 
conservative  local  refinements. 

We  have  developed  a  general  approach  to  maintain  conservation  where  local 
refinements  are  introduced.  The  FAS  coarse  grid  equations  near  the  refinement 
boundary  are  suitably  modified  for  this  purpose:  a  special  term  is  added  to  the 
fine- to- coarse  defect  correction  r. 

Using  such  conservative  refinements  for  the  above  model  problem  gave  the 
accuracy  one  should  expect.  That  is,  the  solution  error  at  distance  r  from  the 
source  is  0(h(r)2/r2),  where  h(r)  is  the  largest  meshsize  at  distances  up  to  r  from 
the  source. 

It  was  also  realized  that  the  general  criteria  for  grid  refinement  [5,  §9.5]  should 
be  modified  in  case  of  interior  singularities,  so  as  to  take  into  account  the  large 
cancellations  among  the  large  contributions  to  the  error  functional  from  the  local 
truncation  errors  'r)  near  the  singularity.  Such  cancellations  are  typically  achieved 
if  conservative  discretization  is  used. 


3.  Multigrid  for  time-dependent  problems 

A  preliminary  description  of  the  results  of  our  contract  work  on  time-dependent 
(primarily  parabolic)  problems  is  contained  in  [  6].  Two  detailed  papers  are  under 
preparation,  containing  results  from  more  extensive  numerical  studies.  It  has  taken 


longer  than  expected,  since  a  student  who  worked  on  this  project  (E.  Gendler) 
could  not  stay  with  us,  and  his  ad-hoc  programs  needed  major  re-writing.  In  fact, 
we  are  still  busy  calculating  results  and  producing  graphs  needed  to  present  a 
comprehensive  and  convincing  account  of  the  new  techniques.  We  will  send  this 
account  to  the  U.S.  Army  as  soon  as  it  is  completed.  Here  we  summarize  the  main 
points;  for  detailed  formulae  and  sample  tests  see  [6]. 


3.1  Multigrid  solvers  for  implicit  systems 


Some  time-dependent  problems  may  need  no  multileveling.  These  are  hyper¬ 
bolic  schemes  where  all  the  characteristic  velocities  are  comparable  to  each  other, 
and  their  explicit  discretization  on  one  grid  is  therefore  fully  effective:  the  amount 
of  processing  is  essentially  equal  to  the  amount  of  physical  information.  However, 
as  soon  as  any  stiffness  enters,  implicit  discretization  and  multigrid  techniques 
become  important. 

The  multigrid  algorithm  can  first  simply  be  used  at  each  time  step  as  a  fast 
solver  for  the  system  of  equations  arising  from  the  implicit  descretization.  This 
use  could  be  considered  quite  straightforward,  but  it  did  produce  some  confusion 
among  several  investigators,  mainly  with  regard  to  the  following  two  questions: 

(i)  In  the  usual  full-multigrid  (FMG)  algorithm  for  steady-state  problems,  the 
first  approximation  to  the  solution  on  a  fine  grid  is  obtained  by  interpolating  from 
a  coarse-grid  approximation  (see  [1,  §4]  or  details  in  [5,  §7]).  In  case  of  evolution 
problems,  should  the  first  approximation  still  be  taken  from  the  coarser  grid,  or 
should  it  simply  be  the  solution  at  the  previous  time  step? 

(ii)  The  FMG  algorithm  in  case  of  steady-state  problems  usually  solves  the 
problem  in  one  cycle;  i.e.,  starting  from  the  above  first  approximation,  one  multi¬ 
grid  correction  cycle  (usually  a  V  cycle;  for  some  problems  a  W  cycle)  is  all 
one  needs  to  obtain  algebraic  errors  small  compared  with  the  truncation  (dis¬ 
cretization)  errors.  Hence  the  algorithm  is  called  1-FMG  (op.cit.).  In  case  of 
time-dependent  problems,  at  each  time  step  one  needs  to  solve  to  the  level  of  the 
incremental  truncation  error,  i.e.,  the  (much  smaller)  discretization  error  added  at 
that  time  step,  not  to  the  level  of  the  accumulated  truncation  error.  Is  the  1-FMG 
powerful  enough  to  do  that? 

The  simultaneous  answer  to  these  two  questions  is  that  the  1-FMG  algorithm 
should  be  consistently  applied  to  the  incremental  problem ,  and  then  it  will  easily 
obtain  algebraic  errors  small  compared  to  the  incremental  truncation  error.  In 
case  of  linear  problems  this  can  be  done  by  writing  the  system  of  equations  as 
equations  for  the  increment  (the  difference  between  the  solution  and  the  solution 
at  the  previous  time  step),  and  applying  the  1-FMG  algorithm  for  this  new  system. 
As  a  result,  for  example,  the  first  approximation  on  the  fine  grid  will  effectively  be 
the  sum  of  the  previous-time  solution  and  an  interpolated-from-the-coarser-grid 


first  approximation  to  the  increment  In  case  of  nonlinear  problems,  the  same  can 
be  done  by  using  the  FAS-type  multigrid  formulation  ([1,  §7]  or  [5,  §8]). 

To  demonstrate  the  above,  we  have  conducted  detailed  tests  with  second- 
order  Crank-Nicolson  discretization  of  the  heat  equations,  with  various  types  of 
initial  conditions:  smooth,  highly  oscillatory,  and  discontinuous.  In  all  cases  a  1- 
FMG  algorithm,  employing  V(l,  1)  (i.e.,  a  V  cycle  with  one  pre- relaxation  and  one 
post-relaxation;  cf.  [1,  §4]  or  [5,  Fig.  1.1])  or  V(2,0)  cycles,  was  used  at  each  time 
step.  We  tested  cases  where  the  errors  can  be  exactly  calculated,  and  showed  that 
the  remaining  algebraic  errors  were  indeed  negligible  compared  to  the  truncation 
errors.  We  have  also  showed  that  some  wrong  multigrid  approaches  produce  fax 
worse  results. 

The  computational  cost  of  the  1-FMG  solver  at  each  time  step  is  several 
(typically  6)  “work  units”,  independently  of  h  and  A:  (the  meshsize  and  the  time- 
step-size),  where  the  work  unit  is  the  computational  work  of  one  explicit  time  step. 
For  many  types  of  problems,  advanced  multigrid  techniques,  discussed  next,  even 
reduce  the  required  work  much  further,  to  a  fraction  of  a  work  unit  per  time  step, 
while  other  multigrid  procedures  (discussed  in  Secs.  3.3  and  3.4)  inexpensively 
enhance  accuracy  and  cure  some  possible  bad  features  of  the  usual  discretization. 


3.2  Frozen  r  techniques 

In  many  problems,  notably  in  parabolic  ones,  at  most  locations  during  most 
of  the  time  the  increment  function  is  very  smooth,  hence  seldom  requires  fine-grid 
processing.  In  the  heat  equation  u<  =  Au  +  F,  for  example,  the  increment  cam 
be  non-smooth  only  near  non-smooth  initial  or  boundary  conditions  or  near  non¬ 
smooth  changes  in  F ,  where  “near”  means  distance  0(h)  in  space  and  0(h2)  in 
time. 

The  frozen- r  technique  is  a  method  to  avoid  unnecessary  fine-grid  processing. 
It  was  described  in  [7,  §3.9].  It  is  based  on  the  FAS  multigrid  formulation  [5 
§8].  The  finest  grid  is  employed  only  once  per  many  times  steps.  When  it  is 
not  employed,  its  r  corrections  to  the  next-coarser-grid  equations  (see  [5,  §8.2]) 
are  calculated  from  previous  times  (assuming  their  decay  in  time  is  proportional 
to  the  decay  shown  in  the  coarser-grid  solution  norm).  The  next  coarser  grid  is 
similarly  skipped:  It  is  employed  only  twice  as  often  as  the  finest;  and  so  on. 
Returning  to  a  grid  after  an  unemployment  interval,  its  values  are  first  calculated 
from  its  former  values  (assuming  again  a  decay  proportional  to  the  decay  in  the 
coarser-grid  solution  norm),  and  then  corrected  from  the  coarser  grid  values  by 
FAS  interpolation  (Eq.  (8.6)  in  [5]). 

Automatic  criteria  to  decide  at  what  time  a  return  should  be  made  to  a  finer 
level  have  been  described  in  [7,  §3.9].  Since  the  criteria  can  be  employed  locally, 
finer  levels  can  be  activated  only  when  and  where  needed. 


In  the  present  research  the  automatic  criteria  were  not  implemented.  Our 
primary  aim  was  to  demonstrate,  through  a  simple  example,  the  huge  reduction 
obtainable  in  the  computational  effort  by  a  uniform  use  of  frozen -r  techniques  for 
uniform  problems.  For  this  purpose  we  chose  the  above  heat  equation  in  a  square 
domain  with  periodic  boundary  conditions  and  with  arbitrary  initial  conditions. 
Since  we  want  to  solve  to  the  level  of  fine-grid  truncation  errors,  the  more  difficult 
case  is  F  =  0,  where  the  solution,  and  hence  also  the  truncation  errors,  tend  to  zero 
as  t  — *  oo.  We  thus  examined  this  case  in  particular.  In  this  case  exact  solutions 
to  the  differential  and  discrete  equations  based  on  Fourier  series  can  be  calculated. 
They  have  been  incorporated  into  our  programs  for  comparison  purposes. 

We  defined  as  our  target  truncation  errors  the  spatial  truncation  error;  i.e.,  the 
difference  between  the  true  differential  solution  and  the  solution  obtained  when  our 
space  discretization  (the  usual  5-point  O(h^)  discretization)  is  used  with  vanishing 
time  step  (the  “method  of  lines”).  This  is  the  error  level  below  which  we  wanted 
to  solve,  by  applying  the  Crank-Nicolson  discretization  together  with  the  frozen  r 
approach. 

Simple  order-of-magnitude  analysis  shows  that,  if  our  finest  meshsize  is  /», 
any  grid  with  meshsize  H  should  be  visited  once  per  time  interval  St  =  0(tH/h 3), 
where  t  is  the  time  from  the  (potentially  non-smooth)  initial  conditions.  Imple¬ 
menting  this  rule  in  our  algorithm,  we  have  indeed  confirmed  that  our  differential 
error  (the  difference  between  our  solution  and  the  differential  solution)  was  always 
comparable  to  the  spatial  truncation  error.  This  has  been  demonstrated  for  both 
smooth  and  non-smooth  initial  conditions. 

From  the  rule  we  implemented,  it  is  readily  seen  that  the  total  work  in  solving 
(to  the  level  of  spatial  truncation  error)  our  evolution  problem  from  t  =  0  to  t  =  T 
requires  only  0(log(T//i2))  work  units. 


3.3  Double  discretization  technique 

The  Crank-Nicolson  (CN)  discretization  is  the  least  expensive  and  most  con¬ 
venient  stable  scheme  with  second  order  accuracy  in  both  the  meshsize  h  and  the 
time  step  k.  It  has,  however,  one  serious  flaw  for  k  /i2:  Transient  highest  spatial 
frequencies  (transient  Fourier  components  with  wavelength  close  to  2 h),  which  in 
the  differential  solution  decay  almost  completely  in  time  k ,  hardly  decay  at  all  at 
each  time  step  (only  their  sign  is  flipped).  In  amother  simple  stable  discretization, 
the  Fully-Implicit  (FI)  scheme,  such  components  properly  decay  in  one  time  step, 
but  accuracy  is  only  first  order  in  k  (hence  k  >  h2  does  not  really  make  sense). 

This  situation  can  be  described  as  a  conflict  between  low  frequencies  (for  which 
CN  is  better)  and  high  ones  (for  which  FI  is  better;  accuracy  order  is  irrelevant 
for  non-smooth  components).  In  multigrid  processing  the  latter  are  treated  only 
at  fine-grid  relaxation,  while  the  former  are  effectively  determined  by  the  coarse- 


grid  corrections.  It  is  therefore  possible  to  resolve  the  conflict  by  using  the  two 
discretization  schemes  at  different  parts  of  the  algorithm  (similar  to  the  “double 
discretization”  schemes  in  steady-state  multigrid  solvers;  see  [5,  §10.2]):  CN  is  used 
in  most  places,  but  it  is  replaced  by  FI  at  relaxation  on  the  finest  grids  during  the 
first  time  steps.  CN  is  still  used,  even  on  each  of  those  finest  grids  even  during  the 
first  time  steps,  for  calculating  the  residuals  transferred  to  the  next  coarser  grid. 

We  have  tested  this  approach  on  our  model  heat  problem,  and  found  the 
expected  behavior:  while  slowly  decaying  high  frequencies  were  completely  elimi¬ 
nated,  the  CN  accuracy  was  always  maintained,  both  for  smooth  and  discontinuous 
initial  data. 

Not  concluded  yet  are  tests  which  combine  this  double  discretization  approach 
with  frozen  r  (Sec.  3.2).  We  expect  this  to  be  a  powerful  combination. 


3.4  r  extrapolation 

In  multigrid  solvers  for  steady-state  problems,  “r  extrapolation”  [5,  §8.4]  is  a 
local  h  extrapolation  based  on  the  fine-to-coarse  correction  rff  appearing  in  the 
FAS  formulation.  It  is  used  to  obtain  a  higher  order  approximation  for  little  extra 
work,  without  employing  higher  order  discretization  formulae.  It  is  more  useful 
than  the  Richardson  extrapolation  since  it  is  local  (extrapolating  the  equation,  not 
the  solution),  hence  it  can  be  used  even  when  global  extrapolation  is  not  available. 
It  can  be  used  for  example  together  with  any  procedure  of  local  refinements. 

We  have  developed  and  tested  a  similar  method  for  time  dependent  problems. 
The  FAS  multigrid  formulation  is  used  as  in  the  frozen  r  technique,  so  its  fine- 
to-coarse  defect  correction  is  calculated,  where  H  =  2h  is  the  coarse-grid 

meshsize  and  K  is  the  time  step  for  which  the  correction  is  calculated.  Choosing 
K  =  2k  in  case  of  CN  discretization,  and  replacing  rh  ’k  by  f7/,  jt  >  can  be  shown 
to  formally  raise  the  approximation  order  from  0(/i2  +  k2)  to  0(/i4  +  fc4). 

We  have  tried  doing  this  for  the  model  heat  problem  described  above.  It  was 
found  that  this  r  extrapolation  indeed  considerably  reduced  truncation  errors,  but 
to  obtain  real  fourth  order  approximation  we  had  to  use  fourth-order  restriction 
(fine-to-coarse  transfer). 

The  reduced  truncation  errors  of  course  makes  the  task  of  solving  to  the  level 
of  these  errors  more  difficult,  hence  requires  modification  in  how  often  a  finer  grid 
should  be  visited  in  a  frozen-r  method  (see  Sec.  3.2).  This  question,  as  well  as  the 
combination  of  r  extrapolations  with  the  CN-FI  double  discretization  technique, 
require  further  investigations. 
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